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PACS 81 . 05 . Rm - Granular Materials 

PACS 61.20.Lc - Time-dependent properties; relaxation 

PACS 64 . 70 . pf - Glass Transitions 

Abstract - We present a mode-coupling theory for the dynamics of a tagged particle in a driven 
granular fluid close to the glass transition. The mean-squared displacement is shown to exhibit 
a plateau indicating structural arrest. In contrast to elastic hard-sphere fluids, which are solely 
controlled by volume fraction, the localisation length as well as the critical dynamics depend on the 
degree of dissipation, parametrized by the coefficient of normal restitution e. Hence the resulting 
glassy structure as well as the critical dynamics are nonuniversal with respect to e. 
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Introduction. — Experimental investigations [IHU as 
well as simulations [3 [6] of dense assemblies of agitated 
granular beads show signature's of a transition from fluid 
to glassy behavior. Motivated by the striking similarity 
etween the measurements for granular systems — neces- 
sarily far from equilibrium — and for equilibrium colloidal 
systems, we have recently generalized the mode-coupling 
theory of the coherent density autocorrelation function 
4> q {t) to driven dissipative systems [7J. Here, we extend 
the analysis to the single-particle dynamics and compute 
the incoherent density autocorrelation function (jf„{t) as 
well as the mean squared displacement (MSD), which are 
directly accessible to experiments. 



. . In two dimensions, experimental measurements of the 
. £h Sr 2 (t) are available by direct imaging from air flu- 

^ idized [5] and mechanically agitated systems fj] . In three 
j_i dimensions, the MSD is observed by diffusive-wave spec- 
.^troscopy (DWS) in gravity-driven flows of glass beads 
in ambient air pQ and in a water- fluidized bed [2]. All 
these studies as well as computer simulations [5J in two- 
dimensional systems reveal the development of a plateau 
in the MSD and relate their findings to glassy dynam- 
ics. With analogous results from colloidal dispersions in 
mind, this has been interpreted as the signature of a gran- 
ular glass transition. Measurements of additional observ- 
ables support this interpretation. Namely, the incoherent 
scattering function has been found to develop a two-step 
relaxation at high densities both in experiments 0] and 




Fig. 1: Incoherent scattering functions 4>q(t) for coefficient 
of restitution e = 1.0 (dashed, elastic case), 0.5 (full curve), 
and 0.0 (dotted) for packing fractions <p from right to left: at 
the glass transition <p c (e), and at 0.999<,o c (£), 0.99ip c (e), and 
0.9v? c (e), respectively. Solutions of Eq. (|la[) are presented for 
the wave vector qd — 4.2. For e = 0.5, an additional solution 
is shown for <p — 0.4 and accompanied with a solution where 
the damping v q is set to zero (chain line). 



in simulations [B]. Moreover, plateaus extending over in- 
creasingly lager windows in time imply a strong decrease 



of the diffusion coefficient D = linn 



Sr 2 (t)/6t G3EJE] 
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accompanied by a strong increase of relaxation times, r, 
4,6 upon approaching the granular glass transition den- 
sity. 

For systems in thermal equilibrium, mode-coupling the- 
ory (MCT) has become an established tool for the investi- 
gation of glassy dynamics, it describes many experimental 
features and has the potential for non-trivial predictions 
[5]. Applied to the case of colloidal suspensions, it is found 
that MCT is quantitatively accurate to about 20% in the 
density. For the mean-squared displacement, MCT de- 
scribes the measured data for the entire regime available 
which is over eight orders of magnitude in time [SlllOj. 

The granular mode coupling theory for homogeneously 
driven systems, outlined in Ref. [7], predicts (1) the exis- 
tence of a glass transition, i.e., lim^oo (j) q {t) = f q > 0, (2) 
an increase of the packing fraction at the transition, y°(e), 
with increased dissipation quantified by the normal coef- 
ficient of restitution e, and (3) changes to the dynamical 
exponents with dissipation. 

While the glass transition within MCT is found to be 
a singularity in the coherent functions 4> q {t), incoherent 
scattering functions (f> s q (t) and in particular the MSD, 
8r 2 (t), have been measured in the experiments and simu- 
lations discussed above. In order to compare to such data, 
the MCT for incoherent functions is derived for granular 
systems in the following. 

Model. — We consider a granular fluid consisting of 
N hard spheres of mass m = 1 and diameter d in a vol- 
ume V. Particle positions and velocities are denoted by 
{fi, Hi} and we will consider the thermodynamic limit such 
that the density n = N/V is finite. Energy dissipation in 
binary collisions is modeled by incomplete normal resti- 
tution, quantified by a constant coefficient of restitution 
e. To achieve a stationary state the system is driven 
randomly and homogeneously: All particles are kicked 
stochastically at random time intervals. This allows the 
system to relax to a stationary state with finite tempera- 
ture, T, defined as the average kinetic energy of the par- 
ticles. 

It is important to note that the system is not in equilib- 
rium so that the TV-particle distribution function, w(T), 
is in general unknown. We assume that positions 
and velocities are uncorrelated, w(T) = w r ({ri})w v ({vi}) , 
and that the velocity distribution factorizes into a product 
of one-particle distributions, w v ({vi}) = Y[i w i(vi). The 
precise from of Wi(v) is not needed, it only has to satisfy 
W Si = ', finite. Furthermore, the system is assumed 
to be isotropic and homogeneous except for the excluded 
volume interaction: w r ({fi}) = Yii<j ^( r y — d). 

Granular Mode-Coupling Theory. The dynam- 
ics can be formulated in terms of a Pseudo-Liouvillc oper- 
ator [TJ] , so that techniques for the derivation of MCT for 
energy conserving systems with Newtonian dynamics in 
|13j can be used with appropriate changes [7j . The central 
quantity to describe the dynamics of a tagged particle is 



the incoherent intermediate scattering function 



n(t)-n(0)) 



where [. . . ] denotes the average over the random driving. 

Using appropriate Mori-projectors, one finds that <fi q (t) 
obeys the following equation of motion 



d 2 ^ s q {t) + v s q d t d> q {t) + 
+Q 2 g f Q dt'm q (t 
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s (t) 



0. 



(la) 



formally identical to the equation of motion for an energy 
conserving system. Here the friction, is given by the 
same expression as for the coherent part: 



1 



v E - 



-[l-jo(qd)+2j 2 (qd)] 



(lb) 



with ji the spherical Bessel functions and with the classi- 
cal (equilibrium) Enskog collision rate ve = 4y/irnd 2 ga^o 
|14j where gd is the contact value of the pair distribution 
function g(r) and vq — y^T/m is the thermal velocity. 
Following Ref. [7], we applied the mode-coupling approx- 
imation to the memory kernel m s q (t), 
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(1c) 

Details of the derivation from the microscopic dynamics 
will be published elsewhere. Here the <f> q {t) are given by 
the solutions of the coherent MCT equations [Jj, S k de- 
notes the static structure factor and c k the direct cor- 
relation function. In contrast to the coherent case, the 
frequency fl 2 q = q 2 T does not carry any dependence on 
the coefficient of restitution e. The equations of motion, 
Eq. dlj, are solved with the initial conditions <fi q (0) = 
l,dt4> q (0) = 0. The numerical algorithms for solving the 
equations of motion for the coherent and the incoherent 
intermediate scattering function as well as for the MSD 
have been introduced previously [Tl fTUllTS] . 

Intermediate Scattering Functions. Fig. [T] 

shows the solutions for the wave vector qd = 4.2 which is 
at about half of the principal peak in the static structure 
factor S q . For increasing packing fraction ip = irnd 3 /6, 
the scattering functions develop a plateau for < e < 1. 
For smaller e, the critical plateau value f q = lim^oo <j) s q {t) 
at the transition point increases. It is seen in Fig. [2] that 
this increase with e applies to all wave vectors for f q . The 
increase results in a growing half width of /* as a function 
of wavenumber, which is related to the inverse localization 
length of a tagged particle. 

One important result of MCT is the existence of two 
timcscales, which both diverge with the distance from the 
critical point, a = (tp c — ip)/ip c . The first one, t a oc 
(j-iA 20 ) ru les the dynamics near the plateau, whereas 
r oc er~ 7 ,7 = l/(2a) + 1/(26) determines the asymptotic 
time dependence, i.e., the a— relaxation and the critical 
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Fig. 2: Glass form factors (upper panel) and critical am- 
plitudes hg (lower panel), cf. Eq. ([2]), for the incoherent scat- 
tering functions for coefficient of restitution e = 1.0 (squares), 
0.5 (circles), and 0.0 (diamonds). The inset shows the expo- 
nent parameter A (*), the von-Schweidler exponent b (+), and 
the packing fraction at the transition ip c ( x ) depending on the 
coefficient of restitution s. 

behavior of the diffusion coefficient. As for the clastic 
case, both exponents a and b are expressed by a single ex- 
ponent parameter A, which varies with e as shown in the 
inset of Fig. g] 

The plateau values, or glass form factors, f q form the 
basis for the asymptotic expansion of the MCT equations 
of motion around the plateau [16] , Close to the transi- 
tion point, the correlation function in the vicinity of the 
plateau can be expanded in leading order as 

ct> s q (t;a) = f° + h s q G <T (t) (2a) 

which defines the critical amplitudes h q shown in the lower 
panel of Fig. [5] While f q and h s are fixed by the details of 
the equations of motion at the transition point, tr = 0, the 
scaling function G a (t) depends only on the time t and the 
distance a from the transition point. The scaling law (12a[) 
reveals a factorization property unique to glassy dynamics, 
where the complex dynamics can be separated into a time- 
dependent and a wavenumber-dependent part. Directly at 
the transition point, Eq. (|2a[) reduces to the critical law 

ct> q {t;a)=r q +h s q {t/to)- a , (2b) 

with a microscopic time scale t$ and the exponent a which 
is shown in the inset of Fig. 2 of [7]. In the fluid state, 



close to the transition, the above power law describes the 
approach to the plateau value. In the fluid state below 
the transition, G a (t) describes the decay from the plateau 
to zero and in this regime is known as the von-Schweidler 
law 

<t> s q {t;a)=r q -h s q {t/T) h . (2c) 

For asymptotically long times, t 3> t a , the incoherent 
scattering functions obey the socalled a-scaling 

^(t;a)=^(t/r(a)) (3) 

which connects to the von-Schweidler law for x = t/r < 1, 
whereas for large x, the decay is close to a Kohlrausch law 
|17j . and crossing over to an exponential for the largest x. 
The scaling suggested by Eq. ([3]) is applied to the results 
for e = 0.5 from Fig. Q] and displayed in Fig. [3] Time 
scales r are determined where the correlation functions 
cross the value <^(t) = 0.4, and the curves are scaled on 
top of the first correlator accordingly. When getting closer 
to the transition, the correlation functions follow a mas- 
ter curve for progressively longer times. In contrast, an 
equally extended correlation function for the elastic case, 
i.e., a different value for e, clearly violates that scaling. 
With the form factors, f q , and amplitudes, h q , as well as 
the exponents b, all being nontrivial functions of e no such 
scaling is expected. 




Fig. 3: a-scaling of the correlators of Fig. [T] for coefficient 
of restitution e = 0.5 (full curves). Curves are matched at 
4>q(t) — 0.4. For comparison, the last fluid curve for e = 1.0 is 
shown dashed. 

The divergence of the time scale, r, determines the van- 
ishing of the diffusion constant D. Both singularities are 
governed by the asymptotic law D a t" 1 cx <t 7 which is 
shown by full lines in Fig. [4] Individual symbols in the 
same plots show the actual values retrieved from the nu- 
merical solution of the equations of motion. The upper 
panel demonstrates that for sufficiently large time scales 
the asymptotic law describes the numerical values satis- 
factorily. On the contrary, for distances around 10% from 
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the critical point and larger, one cannot expect the asymp- 
totic law to hold — and the early part of the divergence 
might even suggest a power law with a different exponent. 



8r 2 (t)/d 2 
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Fig. 5: Mean-squared displacement Sr 2 (t) at the glass- 
transition packing fraction tp c (e) for coefficient of restitution 
e — f.0, 0.5, and 0.0 decreasing from top to bottom at long 
times. Full curves show the numerical solutions of Eq. @. 
The dashed curve displays the critical law in Eq. (|5c|) with 
to — 0.035 for the case e = f .0. The inset shows the localiza- 
tion length, cf. Eq. (|5a|) . as a function of e with e = f .0, 0.5, 
and 0.0 shown as filled circles. The chain curve shows the elas- 
tic limit r c (e — 1) scaled by the evolution of the mean particle 
separation, £o oc 1/ ^/tp c (e)/tp c (l). The dashed curve indicates 
the calculated localization length with the latter dependence 
scaled out. 



l — q Sr (t) + 0(q ) for small wavenumbers. The equation 
of motion for the MSD in the granular case reads 



Fig. 4: Dependence on the packing fraction, cp, of the time 

scale r where the incoherent scattering function 4>q(r) = 0.1 and can also be obtained from the expansion of (j>g(t) 
(upper panel), and diffusion coefficients D (lower panel) de- 
rived from the numerical solutions of the equations of motion 
for coefficient of restitution e = f.0 (squares), 0.5 (circles), and 
0.0 (diamonds). The full lines show the corresponding asymp- 
totic laws r oc [ip c (e) — ip(e)]~~' and D oc [tfi c (e) — co(e)] 7 . For 
e = f.0, 0.5, and 0.0, the values for 7 are 2.46239, 2.3492f , and 
2.28282, respectively. 



d t 5r 2 {t) + ^u E 8r 2 {t) 

+v 2 f* dt'm^{t~t')5r 2 {t') 



6v 2 t 



(4a) 



The evolution of the diffusion coefficients in the lower 
panel of Fig. @] suggests a quite similar behavior as for 
t. Again the asymptotic law works well very close to 
the transition while for larger distances of 10% and more, 
the numerical solution diverges markedly slower than the 
asymptotic law would suggest. 

The structure of the Equations (|2I3[) is the same as for 
the incoherent functions of the elastic hard-sphere system 
[16| , including the factorization property for the dynamics 
around the plateau. However, the exponent parameter A 
as well as glass form factor (/|, cf. Fig. [2]) and the critical 
amplitude (h q , cf. Fig [2]) do depend on the coefficient of 
restitution. Hence we conclude that the dynamics is not 
universal with respect to dissipation. 

Mean- Squared Displacement. Experiments as 
well as simulations focus on the MSD which is defined 
by 



Sr 2 (t) :-- 



r!(0)) 2 



and the memory kernel within MCT is given by m(°)(t) = 
limg-^o q 2 m q (t), and reads 



1 



e n 



2 6tt 2 



dk k A S k c 2 k MtWk(t) ■ (4b) 



The behavior of the MSD at the respective critical 
points for several e is shown in Fig. [5] The effective En- 
skog damping coefficient 1^(1 + e)/2 decreases for smaller 
e, and this lesser damping causes the MSD for e = 1 to be 
slightly larger than that for e = 0.5 and e = 1 in the time 
window —2 < log 10 t < —1 after the ballistic regime. For 
macroscopic times the solutions reverse their order and fi- 
nally reach their long-time limits 6r 2 with the localization 
length being defined by 



r 2 c = lim 5r 2 (t)/6 



1/ lim m 



(0) 



(*)• 



(5a) 



Similarly to the glass-form factors, for the MSD an asymp- 
totic expansion can be performed with the result similar 
to Eq. © 

\t)/6 = ri-h M suG <T (t). (5b) 
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Together with the critical, 

5r 2 (t)/6 = r 2 c -h MSD (t/t )- a , (5c) 

and the von-Schwcidlcr law, 

6r 2 (t)/6 = r 2 + h MSU (t/T) b , (5d) 

the analysis of the elastic hard-sphere systems in [16] car- 
ries over to the dissipative case. The evolution of r c and 
^msd with e are shown in the insets of Figs. [5] and [5] In 
addition, the limited applicability of the asymptotic criti- 
cal law due to large corrections to scaling applies to elastic 
as well as dissipative hard spheres [TO] . 
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Fig. 6: Mean-squared displacement Sr 2 (t) for coefficient of 
restitution e = 0.5 and packing fractions tp = ip c , 1.001<p c , 
1.01v? c , l.l<p c , on the arrested side — 0.999<p c , 0.99ip c , 0.9p c , 
and 0.4, respectively, on the diffusive side. The terms liquid 
and glass indicate the diffusive and arrested regimes, respec- 
tively. Dashed curves show the results for the critical point 
and tp = 0.4 with the choice of v q = v s q = ve = 0. The initial 
ballistic law v$t 2 is shown by a dotted straight line. The inset 
exhibits the critical amplitude /imsd as a function of s. 

It is seen in the inset of Fig. \b\ that r c decreases with 
e which is corresponding to a higher packing fraction at 
the glass transition (cf. inset of Fig. [5]). While the overall 
decrease of r c from e = 1 to e = is more than 16%, 
the trivial contribution from the increase in density can 
explain only 4% of the decrease: A length scale can be 
defined by the inverse cubic root of the packing fraction 
which gives a scaling factor y/ tp c (e) / tp c (e = 1). When the 
calculated results are scaled with this factor, the dashed 
curve in the inset of Fig.[5]is found. Hence, the major part 
of the decrease of the localization length is a more involved 
prediction than just a simple density scaling argument. 

Short-Time Dynamics. — It is known that the 
damping v q together with the early part of the memory 
kernel can overestimate the total damping considerably, 
which is an issue that cannot be fully resolved [5]. It is 



possible to estimate the order of magnitude of that effect 
without introducing mathematical inconsistencies by set- 
ting all the damping terms, v q , v q , v& to zero and solve 
the equations. The result is shown as the chain curve 
for 4> q (t) in Fig. [TJfor e — 0.5 and tp = 0.4. Outside the 
transient regime, the undamped solution is shifted by 45% 
compared to the damped solution. For the MSD in Fig. [5] 
for the same values of e and tp, it results in a similar shift 
of 45% in time scales. At the transition point, however, 
the solutions with and without damping deviate from each 
other for the relatively large window —2 < log 10 t < 1. It 
is also seen that the curves with Enskog damping devi- 
ate from the short-time asymptote Sv 2 t 2 considerably ear- 
lier in time. Together, those effects can mask the glassy 
dynamics expected for moderately large windows in time 
that are accessible to experiments and simulation. 

Summary and Outlook. — In conclusion, the results 
discussed above and in [7] suggest for testing the theory 
experimentally. While for the variation of the packing 
fraction, tp, one expects to recover most features known 
from the thermal glass transition — from the variation of 
£ the following scenario should emerge: 

(1) The glass transition shifts to higher packing frac- 
tions tp c for lower e, cf. inset of Fig. [5] Since the overall 
change is around 10% it should be measurable directly. 
The behavior of time scales and diffusion coefficient of the 
MSD demonstrated in Fig. [4] supports the possibility of 
extracting this change from a dynamical experiment; even 
the asymptotic scaling law can be used reliably. For a dis- 
cussion of experimental data in three dimensions one must 
first consider that the value for the hard sphere system in 
its elastic limit can be assumed at 92° sa 0.58 [8lll"0]. 
Hence, also the transition at smaller e will be shifted to 
yet higher values of tp. 

(2) The localization length r c decreases for the mean- 
squared displacement (cf. inset of Fig. [5]), and the corre- 
sponding plateau values /* (cf. Fig. [5]) increase for smaller 
£. While an absolute determination of f q and r c from mea- 
sured Sr 2 (t) and <fi q (t) may be ambiguous due to limited 
windows in time, a comparative measurement between dif- 
ferent £ should be within experimental resolution for the 
MSD, cf. Fig. El 

(3) The a-scaling function for the long-time decay of the 
correlation functions and the MSD changes with e. Since 
absolute changes in the critical exponents are comparably 
small, see inset of Fig. [5J direct observation of changes 
in the exponents is probably rather difficult. However, a 
qualitative method is given by testing the a-scaling: For 
fixed £, curves for different tp scale along their plateau 
value onto a single master curve, cf. Fig. [TJ and [3] In 
contrast, when going on a path along the transition in the 
£ direction, a violation of this a-scaling is expected. 

While one can expect that the overall behavior of the 
MCT for granular systems is similar in two and three di- 
mensions as is the case for the elastic hard-sphere system 
[T5] . the necessary glass form factors, critical amplitudes, 
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and exponents cannot be estimated without performing [19] Durian D. J., J. Phys.: Condens. Matter, 12 (2000) 
the actual calculation. The dynamical scenarios observed A507. 
in two dimensions [21E] nevertheless support the existence 
of a glass transition in a driven granular fluid. In addi- 
tion, a scenario for the evolution of <j> s q (t) is found in [3] 
that is quite reminiscent of the 3D case shown in Fig. [TJ 
However, in contrast to the binary mixture in [3J, the re- 
sults from [3] have to be taken with some caution as the 
monodisperse experiment has a higher tendency towards 
ordering which may influence the observed dynamics. In 
all two-dimensional experiments, the estimated plateau of 
the MSD is around 4r^ w O.OlcP, which is consistent with 
results anticipated from Fig. [5] However, the values re- 
ported for the localization length in Refs. pQJJ] are smaller 
than predicted here by a factor of 10 3 . One possible expla- 
nation for such a deviation is the influence of accelerations 
due to gravity |19j which might be remedied by future ex- 
periments under microgravity. 

* * * 
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